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Abstract 



We study hadron properties near the deconfining transition in the quenched lattice 
QCD simulation. This paper focuses on the heavy quarkonium states, such as J/tp meson. 
In order to treat heavy quarks at T > 0, we adopt the 0(a) improved Wilson action on 
anisotropic lattice. We discuss cc bound state observing the wave function and compare 
the meson correlators at above and below T c . Although we find a large change of correlator 
near the T c , the strong spatial correlation which is almost the same as confinement phase 
Oh, survives even T ~ 1.5T C . 



1 Introduction 



It is generally believed that the Quantum Choromodynamics (QCD) exhibits a phase transi- 
tion at some temperature T c , and quarks and gluons confined in the low temperature phase 
are liberated to form the "quark gluon plasma". At the beginning of 2000, CERN reported 
that the QGP state had been created in the heavy ion collision experiment^. In these exper- 
iments, the J/tj) suppression M is regarded as the key signal of QGP formation. Since c-quark 
is heavy ( m c = 1.15 ~ 1.35 GeV||) cc pair are hard to be generated except by the primary 
collisions of nucleons in the high energy heavy ion, i.e. the J/i/j is not created with the thermal 
effects after the QGP formation. Therefor it is expected that the effects of deconfining clearly 
appear in this signal. Further investigation will be performed in RHIC project at BNL. 

On the theoretical side, in spite of the various approaches |fj, |5], q], we are still far from the 
definite understanding of hadron properties near the transition and the fate of the hadronic 
states in the plasma phase. Since the phase transition changes the relevant degree of freedom 
of the system, the model approach which a priori assumes dynamical degrees of the system 
is difficult to treat physics near the phase transition. We investigate these problems using 
lattice QCD, which enables us to incorporate the nonperturbative effect of QCD from the first 
principle. 
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For a long time, in the lattice QCD simulations, hadron masses at finite temperature have 
been argued with the spatial correlation (screening mass) . On the other hand, the study of 
the temporal correlation, which is related to the pole mass, has been started rather recently |8|. 

The previous work|8| studied hadron properties near the T c based on light quarks. They 
caught a sign of chiral symmetry restoration above T c , and a change of correlators in temporal 
and spatial direction near the T c . From the discussion of wave function, however, they find the 
same strong spatial correlation at T > T c as that of below T c , in other words hadronic mode 
at a long range survives in the deconfinement phase. At present, there is no well-established 
way of lattice simulations to attack the spectroscopy at T > 0. However, our recent work 
seems one of the best approaches to the problem. It is interesting to apply this analysis to 
the heavy quarks. In this paper we focus on the heavy quarkonium state, which plays an 
important role as a signal of the quark gluon plasma formation @, B|. 

Our goal of the investigation is the prediction about J ftp suppression or mass shift of 
charmonium Q as signals of QGP formation and the understanding of hadron properties and 
nature of QGP phase near the transition. 

For the study of charmonium physics at T > on a lattice, there are several problems. 
Then we classify these problems into two category, and discuss them individually. 

(i) Precise calculation of temporal correlator of charmonium at T > 

(ii) Extraction of physical properties of charmonium from the correlators 

Firstly we consider the former one. In the lattice QCD simulation at T > 0, we set a temporal 
lattice extent to 1/T. At high temperature, one needs the large lattice cutoff to work with 
the sufficient degrees of freedom in the temporal direction. In order to obtain the detailed 
information of temporal meson correlators at T > 0, a high resolution in temporal direction is 
needed. The large lattice cutoff is also necessary to study a correlator of meson with the heavy 
quarks because of its rapid decreasing behavior. If one tries to overcome these difficulties with 
straightforward way, the tremendous large computational power is necessary. In order to get 
the sufficiently fine resolution with limited computer resources, we use the anisotropic lattice, 
which has a finer temporal lattice spacing a T than the spatial one a a . 

In this work we adopt the same strategy as the previous work|8[ which was tractable to 
analyze light hadrons at T > 0. This strategy is as follows. Firstly the mesonic operator 
are defined, then we observe its correlator and wave function (in the Coulomb gauge) at 
T = 0. Next we investigate how they are affected by the temperature. In order to investigate 
the temperature effects for the state of interest, for example a ground and excited state of 
charmonium, we have to make the good mesonic operator which has the large overlap with 
the state. Because we should compare correlators at T > within shorter temporal lattice 
extent. The wave function gives hints for existence of mesonic state at T > 0. Especially we 
are interested in that of deconfinement phase. In the case of charmonium this wave function 
is an important quantity concerning the J ftp suppression. 

We prepare two sets of gauge configuration whose lattice spacing are different. We con- 
trol the temperature by changing temporal lattice extent Nf. Then above investigation is 
performed on these configurations. 

This paper is organized as follows. In the next section we define the quark action on 
anisotropic lattice and discuss the dispersion relation of free quark and calibration of quark 
field. Sect. |3| is the preparation for study of charmonium correlator. Here various parameters 
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of gauge configurations are determined and the calibration are performed. In Sect. ^| we 
report the charmonium spectroscopy and construction of optimized operator using variational 
analysis at T = 0. Sect. || describes correlators at T > and the measurement of the wave 
function and compare these with results at T = 0. The last section is the conclusion and 
discussion. 



2 Quark action on the anisotropic lattice 
2.1 Quark action 

To treat the quark field on the lattice, we adopt the 0(a) improved Wilson quark formulation. 
To construct the quark action on the anisotropic lattice, we follow El-Khadra, Kronfeld and 
Mackenzie [|l0||, for the following advantages. They expand the lattice Hamiltonian in the power 
of a, and determine the coefficient of each operator by matching the lattice Hamiltonian with 
the the continuum one except the redundant operators. The resultant action takes the same 
form as the clover quark action [pAj in the limit of m — ► 0. On the other hand, in the heavy 
quark mass region (m q 3> Aqcd), the effective-theoretical treatment of the quark action 
enables us to use it for such a quark on a lattice of moderate cutoff. Since the full quark mass 
dependence is incorporated, the same form also covers the intermediate quark mass region, 
and then the small and the large mass regions are smoothly connected. Although our main 
target in this paper is the charm quark, it would be useful to take the other mass region 
into account for the future applications. In addition to these advantages, their argument 
is naturally in accord with the anisotropic lattice. They introduced the different hopping 
parameters for the spatial and the temporal directions, and proposed to tune them so that 
the rest and the kinetic mass take the same value. Such a treatment is inevitably required on 
the anisotropic lattice to assure that the anisotropy of the quark and the gauge fields coincide, 
especially if one employ the dispersion relation for the definition of the quark field anisotropy. 



On the anisotropic lattice, the quark action takes almost the same form as in the Ref. [1C] 

S F = ^(x)K(x,y)^(y), (1) 

x,y 

K(x,y) = 5 Xty -K T [(l- l4 )U 4 (x)5 x+ly + (l + -f 4 )ul(x-4)5 x _ ly } 
-k* E {( r " vM x ) S x +i,y + ( r + 7i)Uj(x - i)6 x _ ly } 

i 

rK a c B E (x)Sx,y (2) 

The spatial and the temporal hopping parameters, K a and k t respectively, are related to the 
bare quark mass itlq and the bare anisotropy parameter jp as follows. 

Ka = 2(m +7F + 3r)' (3) 

K T = ~f F K a . (4) 

where the mo is in the spatial lattice unit. In the free quark case, the bare anisotropy 
is taken to be the same value as the cutoff anisotropy £ = a a /a T . In practical simulation, 

1 The notation in this paper is slightly different from the Ref. |l0[ . 
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the anisotropy parameter receives the quantum effect and should be tuned to give the same 
renormalized anisotropy for the fermion and the gauge fields. This "calibration" will be 
described later. 

There have been used two choices of the value of the Wilson parameter r for the anisotropic 
0(a) improved quark action. In this work, we adopt the choice^] r = l/£. In this case, the 
temporal and the spatial directions are treated in the equal manner in the physical unit. 
As the result, the tree level dispersion relation holds the axis-interchange symmetry in the 
lowest order of p 2 . On the other hand, this choice decreases the masses of doublers which are 
introduced by the Wilson term to eliminate the unwanted poles at the edges of the Brillouin 
zone. The dispersion relation is examined in the later part of this section. Alternative choice, 
r = 1, are adopted in the Ref. |lf, |J, 0]. In this case, the contribution of doublers would 



not cause any problem, in the cost of manifest axis-interchange symmetry. Since we aim 
to develop the form applicable to the whole quark mass region, r = seems preferable 

especially in the light quark mass region. 

Here it is useful to define k so that which has the same relation with the bare quark mass 
mo as the isotropic case: 

- = — -2( 7F + 3r-4) (=2(m + 4)). (5) 

For the light quark systems, the extrapolation to the chiral limit would be performed in 1/k. 
The coefficients of the clover terms, ce and cb, depend on the Wilson parameter r. In our 
choice r = eg and eg are unity at the tree level. 

We apply the mean-field improvement proposed in the Ref. |l5| . On the anisotropic lattice, 
the mean-filed values of the spatial link variable u a and the temporal one u T are different from 
each other. The improvement is achieved by rescaling the link variable as Ui(x) — ► Ut(x)/u a 
and U^(x) — > U4(x)/u T . This replacement leads the following values for the coefficient of the 
clover terms. ^ 

ce = 5-, cb = -o-. (6) 

The determination of mean-field values of the link variable u a and u T are described in the 
next section. 

In this paper, our target mass region is around the charm quark mass. The temporal 
cutoffs in this work are 4.5 and 6.4 GeV, and well above the charm quark mass. For these 
quark mass and a" 1 , the effective-theoretical treatment would not be necessary to be applied. 
Such consideration will be called for the calculations containing the 6-quark on the same size 
of lattice. In the effective-theoretical treatment, the ratio of the spatial and the temporal 
hopping parameter is tuned so that they give correct dispersion relation of the nonrelativistic 
quark flcfl . On the anisotropic lattice, the calibration automatically incorporates this condition 
if one use the nonrelativistic dispersion relation as the anisotropy condition. 



2.2 Dispersion relation of free quark 

Now we consider how the dispersion relation of the free quark is changed by the introduction 
of anisotropy. Observing the action (|2[), one notices that the larger anisotropy £ causes 
the smaller spatial Wilson term. Then the question is how the contribution of the doubler 
eliminated by the Wilson term becomes significant. The action (g) leads the free quark 
propagator, 

S(P) = — —2, (7) 

Z74 sinp4 + i £7 • S + m t + (1 — COSP4) + 2 r (fi 
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Figure 1: Dispersion relation of anisotropic free quark. Left figure is of Set-I, right one is of 
Set-II. 



where C = £ > m t = w-oC; Si = sinpj and pi = 2 sin(pj/2). Then the dispersion relation of 
the free quark is 

cosh £ ( a = 1+ ^ + ( "" + K' 2 . (8) 
2(1 + m t + \rQp ) 

Neglecting the higher order terms in p and in a, the relativistic dispersion relation E 2 = 
m 2 + C^p 2 holds for the small quark mass, (m 2 = m 2 (l — mt), and mt is the bare parameter.) 

Fig. [j] shows the dispersion relation (Q) at £ = 5.3 and 4.0 for various values of mt. Now let 
us consider the practical cases that a~ l = 4.5 GeV for £ = 5.3 (Set-I) and a~ l = 6.4 GeV for 
£ = 4.0 (Set-II). These values are obtained in our numerical simulations, and described in the 
next section. In the heavy quarkonium, the typical energy and momentum exchanged inside 
the meson are in the order of mv 2 and mv respectively[16|. For the charmonium, v 2 ~ 0.3, 
then typical scale of the kinetic energy is around 500 MeV. It is noted that mt ~ 0.3 for Set-I 
and mt ~ 0.2 for Set-II correspond to the charm quark mass. Let us consider two quarks 
inside meson with opposite momenta p = ±ir/a. Then 2(E(p z = a/ir) — E(0)) ~ 0.5 GeV and 
~ 1 GeV for Set-I and Set-II respectively. Although Set-I lattice may not be free from the 
systematic effect, Set-II would be successfully applicable to the low-lying charmonium system. 

For comparison, we also examine the light quark mass region. For Set-I, m = 0.02-0.06 
corresponds to 90-270 MeV, which is used in the Ref. H as the light quark mass region 
with the anisotropic Wilson quark action with r = 1. E(p) — E(0) rapidly decrease at the 
edge of the Brillouin zone, and the height at z = a/n is around 300 MeV. For two quarks 
with momenta p = ±a/n, additional energy of doublers is ~ 600 MeV, and again seems not 
sufficiently large compared with the typical energy scale transfered inside mesons. In the case 
of Set-II, this value increases to 1.4 GeV, and seems to be applicable to the meson systems 
without large systematic effect. 

2.3 Calibration of quark field 

On anisotropic lattice, the anisotropy of quark configuration £p must be equal to that of gauge 
field f . 

£(7G, If) = Cf(7g, If) (9) 
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Since £ and £p are function of 7g and 7^ in general, a nonperturbative determination of 
the combination of 7g and jp which satisfy the condition @ requires much effort. In the 
quenched case, however, these determination are rather easy to be performed, because £ can 
be determined independently of jp. After the determination of £, one can tune 7^ so that 
the certain observable satisfies the condition @ . We call this procedure as "calibration" . 

There are several determinations of £p. In the Ref. []|], the ratio of the temporal to 
the spatial meson masses, £p = m a /m T , is used as such a observable. However this is not 
suitable for the present case, since the charm quark mass is larger than or comparable with 
the spatial lattice cutoffs. In order to determine £f, we use the dispersion relation of the free 
meson [EjL |i~3| , |l4| . For a heavy quark, one may use the nonrelativistic dispersion relation, 
E = m + p 2 /2m^p. In this paper, we alternatively use the relativistic dispersion relation of 
meson for the calibration. This form is also available for the light quark mass region. 

We assume that the meson is described by the following lattice Klein-Gordon action. 

S = E ^r4> ] {x) \-tlDt - D 2 + m 2 } 4>{x), (10) 
where mo is in the unit of a a . Then the free meson satisfies the dispersion relation 



1 ,-2 

?2 



coshE(p) = \ + --^{f + m 2 ). (11) 



Using this relation, one can determine the anisotropy £p as 



->2 



^ F 2(cosh£(£)-cosh.E(0)) ^ 



This condition forces the axis-interchange symmetry to the meson field. 

3 Lattice setup 

3.1 Gauge configuration 



The numerical calculations performed on the two sets of anisotropic lattices [17], one is with 
the standard plaquette action and the other is with Symanzik type improved action at the 
tree level [|18|. These actions are represented with the following form; 

S G = — E fai(l-Plj(x)) + (*2(2-Rij(x)-R i i(x))] 

^ G x,i<j<3 

+PlG E " P ^ X )) + C 12( 2 " Bui*) ~ R ±ii x ))} , ( 13 ) 

x,i<3 

where the plaquette P fll/ (x) and rectangular loop R fiu (x) are defined as follows, 



P, 



\rbTi[U^{x)U v {x + p)uUx + u)Ut{x)} (14) 



Rp, = -ReTrp^U^x + p)U v (x + 2£)tfJ(x + u + £)£/t(x + v)Ul(x)}. (15) 

The standard action is the case with c\\ = 1 and c\2 = 0, and the improved action is 
en = 5/3 and c±2 = —1/12. The Set-I is the same configurations as used in the Ref. ||. 
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Set 


en 


Cl2 


size (3 


1G 


# conf. 


Set-I 


1 





12 3 x N t 5.68 


4.00 


60 


Set-II 


5/3 


-1/12 


16 2 x 24 x iV t 4.56 


3.45 


120 



Table 1: Simulation parameters for the gauge configurations. 



Set 




a-\GeV) 


a" 1 (GeV) 




N t (T/T c ) 




Set-I 
Set-II 


5.3(1) 
3.950(18) 


0.85(3) 
1.610(14) 


4.5(2) 
6.359(62) 


72 (~ 0), 
96 (~ 0), 


20 (0.93), 16 (1.15), 
26 (0.93), 22 (1.10), 


12 (1.5) 
16 (1.52) 



Table 2: Scale parameters for the gauge configuration. The £p, a~ and a~ are measured 
on the configurations at T = ( Nt = 72 (Set-I), 96 (Set-II)). Here the temperatures are 
roughly estimated in the unit of T c . 

The parameters for these configurations are summarized in Table |l|. These parameters are 
adopted so that the spatial and temporal lattice extent are about 3 fm respectively. 

Both numerical calculations are done on the quenched configurations which are generated 
by the pseudo-heat bath algorithm with 20,000 thermalization sweeps, the configurations 
being separated by 2,000 sweeps. These configurations are fixed to Coulomb gauge. The 
statistical errors are estimated using the jackknife method unless mentioned explicitly. 

We determine the parameters for the gauge field configurations, renormalized anisotropy 
£(= a a /a T ) and spatial cutoff scale a" 1 . We define the ratios of the Wilson loop on spatial- 
spatial (era) and spatial-temporal (err) plane flT9|, |20| , pT| | . 

R a (r,x) = W aa (r,x)/W aa (r + l,x), (16) 
R T (r,t) = W aT (r,t)/W aT (r + l,t) (17) 

Then £ is determined so that the matching condition, R a (r,x) = R T (r,t = is satisfied. 
From this analysis we conclude £ = 5.3(1) (Set-I) and £ = 3.950(18) (Set-II). Cutoff scales 
are determined from the static quark potential using the physical value of string tension 
^/o^I = 427MeV|§2|. The spatial cutoff scales are a" 1 = 0.85(3) GeV (Set-I) and 1.610(14) 
GeV (Set-II) respectively. 

From the calculation of Polyakov loops and its susceptibilities, we find that the tempera- 
tures of lattices with Nt = 18 (Set-I) and 24 (Set-II) are both just above T c . Then we estimate 
T c = 250 ~ 260 MeV from the temporal cutoff scale on our lattices. These estimation are 
consistent with the other quenched lattice results. We prepare the finite temperature gauge 
configurations at just below T c , just above T c and about 1.5T C for each sets. These results are 
summarized in Table [2|. 

Here we notice that the hadronic correlator depends on which the sector Polyakov loop 
stays in Z3 space at T > T c . Since we treat the quenched QCD as an approximation of the 
full QCD, we choose the real sector for the value of Polyakov loop at T > T c . 
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Set 




u T 


CE 


CB 


Set-I 


0.75050(16) 


0.992436(13) 


1.7889 


2.3656 


Set-II 


0.812354(92) 


0.9900788(90) 


1.5305 


1.8654 



Table 3: Mean-field value of link variable and mean-field improved clover coefficients for each 
sets. 

3.2 Mean-field value on the anisotropic lattice 

We state the calculation of mean-filed value on the anisotropic lattice for the mean-field 
improvement of the clover coefficients. 

There are two commonly used methods. One is determined from the expectation value of 
plaquette. This is widely used for its easiness to measure. On anisotropic lattice these are 
determined as follows. 

u a = -{ReTrP ij (x))*, u T = ^(ReTrP i4 (x))5 (18) 
The other is the trace of link calculated in the Landau gauge. 

u a = -{ReTrUi(x)), u T = -{ReTrU 4 (x)) (19) 

In the former case u T is greater than the unity on our lattice. Then the latter definition seems 
more reasonable. In this work we determined the mean-field values in the Landau gauge. 
The Landau gauge fixing is realized by maximization of Eq. ( |20| ) . 

^ReTr {U^x) + U 2 (x) + U 3 (x) + c st U 4 (x)) (20) 

X 

Here, the temporal coefficient c s t appears on the anisotropic lattice. We adopt the self- 
consistent mean-field improvement of c s t ■ As the tree level we chose c s t = £, 2 ■ Using the mean- 
field improved c s t = ^u a ju T we calculate the mean-field values recursively. We perform these 
calculation with 20 configurations. In our case, the result of 3rd measurement is consistent 
with the input u a /u T which is determined from the linear interpolation of the tree level and 
first mean-field improvement result. These results are summarized in Table || together with 
the mean-field improved ce and cb- Here we mention that the mean-field improved 7g give 
the reasonable estimation for £ within 1% (Set-I) and 6% (Set-II) error. 

3.3 Calibration result 

We now turn to the calibration of quark field. It is performed by the Eq.(|i~l|) with the mo- 
mentum p = (0, 0, 0) and (0, 0, 1) for each mesonic channel, pseudoscalar (Ps) and vector (V). 
These results are summarized in Fig. |2|. These calibrations are performed at the parameters 
which correspond to the pseudoscalar mass mp s ~ 2.6, 3.0 and 3.4 GeV respectively. 

The detailed definition of meson correlator which is adopted here is shown in Sect. 4.1 
In these calculations we use the smeared source function which is determined from the wave 
function at 1/k = 9.868 (Set-I) and k = 8.797 (Set-II). Here we mention that the shape of 
wave function has very weak k dependence. 
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3.6 3.8 4.0 3.6 3.8 4.0 3.6 3.8 4.0 ' 3.1 3.2 3.3 3.1 3.2 3.3 3.1 3.2 3.3 

Yf Yf Yf Yf Yf Yf 



Figure 2: The left figure shows the calibration for Set-I and the right one is for Set-II. In the 
both figures, values of £ ( =5.3(1) (Set-I) and =3.950(18) (Set-II)) are indicated. The three 
values of k correspond to a pseudoscalar mass m Ps ~ 2.6, 3.0 and 3.4 GeV respectively. 





1/k 




IF 




10.300 


0.093545 


3.629(91) 


Set-I 


9.868 


0.096098 


3.703(84) 




9.480 


0.098599 


3.765(79) 




9.041 


0.110350 


3.251(52) 


Set-II 


8.797 


0.113121 


3.262(49) 




8.590 


0.115564 


3.272(47) 



Table 4: Calibration results for Set-I and Set-II. These parameters satisfy the condition of 
£ = £p within the statistical error for the both mesonic channel. The error of 7_p is estimated 
from the error of £ and £p 

From the Fig. |, we find that the 7^ dependence of £p is almost linear and k dependence 
of its slope is small. The 7^ which is determined from the pseudoscalar and the vector gives 
consistent value within the statistical error for each k. Especially £p for Set-II are in good 
agreement with the pseudoscalar and the vector. We determine the 7^ which satisfies the 
condition of Eq.@ from the interpolation with the linear fitting. The error of 7^ is estimated 
from the error of £ and £p. These calibration results are summarized in Table ||. Here we 
mention that the mean-field improved 7^ give the reasonable estimation for £p within 9% 
(Set-I) and 1% (Set-II) error. 

4 Results at Zero temperature 
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4.1 Meson correlators 



The mesonic operator which is used in this paper is in the following form, 

0^\z,t) = J2u(y)q&t) lM q(z + y,t), (21) 
y 

where u>(y) is a smearing function and 7 m = 75 and 71 for the pseudoscalar and the vector 
respectively. Using these mesonic operators we construct the meson correlators G M ^ w (t) in 
the Coulomb gauge as the following, 

Gj&^W = £<0if (0)) (22) 

z 

According to our strategy, we need to construct a good operator which has large overlap with 
the ground state ( or the excited state ) of heavy quarkonium. For this purpose we examine 
two types of correlators. 

In the first type of correlator we use the smeared source and point sink ( u/(y) oc 5{y) 
), which is already used in the calibration described in Sect. 3.3 . This source function is 
defined from the measured wave function u;(j?)||23[| so that the smearing reflects the actual 
distribution of quark and anti-quark. This smearing procedure works well for the suppressing 
higher excited state to the correlator. The wave function is well fitted to the form 

w(y) oc exp (—a\y\ p ). (23) 

The determined parameters a and p are used to generate the smearing function. These 
correlators are used in the charmonium spectroscopy in the next subsection. 

For the study at T > we need more systematic optimization of the correlators. We apply 
the variational analysis, and regard the diagonalized correlator as the optimized correlator. 
This analysis can extract not only the ground state but also the excited state. In this analysis 
we use the correlators with smeared source and sink, whose smearing functions are determined 
by solving the Schrodinger equation with the potential model. This technique is mentioned 
in Sect. 4.3 . 

4.2 Charmonium spectroscopy 

We show the results of spectroscopy for the heavy quarkonium at T = 0, which is the basis of 
the study at T > 0. In this section we also examine how our action works in the heavy quark 
system. The pseudoscalar and the vector meson channels are measured with the parameters 
listed in Table |. Fi gure y show the effective mass plot at T = 0. In this calculation we 
use the correlators with the smeared source introduced in the previous subsection. Here the 
parameters for the smearing function are a = 0.7096, p = 1.1267 for Set-I and a = 0.2978, p = 
1.1918 for Set-II. The plateaus of the effective mass appear beyond t ~ 10 (Set-I), 20 (Set- 
II) . Concerning the determination of smearing function our choice works well in heavy quark 
system. On the other hand in the light quark system the wider smearing function enhances 
the ground state contribution ||. In the small t region, contributions from the excited states 
remain. However this small t region is the main stage of the study at T > 0. In order to 
estimate the temperature effects of the ground state contribution at the high temperature, we 
need further improvement of the mesonic operator. In the next subsection more systematic 
study with variational analysis is examined. 
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Figure 3: Effective masses at T = 0. The left figure is for Set-I and the right one is Set-II. 



1/k 


m Ps 


m v 


HFS 


10.300 
Set-I 9.868 
9.480 


0.77934(75) 
0.68235(77) 
0.59381(78) 


0.79146(83) 
0.69600(88) 
0.60917(91) 


0.01213(26) 
0.01365(30) 
0.01535(35) 


9.041 
Set-II 8.797 
8.590 


0.54170(39) 
0.47236(39) 
0.41029(40) 


0.55050(43) 
0.48238(44) 
0.42182(45) 


0.00880(15) 
0.01002(16) 
0.01152(18) 



Table 5: Spectroscopy for Set-I and Set-II with parameters in Table |I[ HFS means the 
hyperfine splitting for a heavy quarkonium, m Ps — m v . There results are shown in a T unit. 



We fit the correlator to the single exponential form at t = 30-36 (Set-I), 40-48 (Set-II). 
These results are summarized in Table [f| The parameters with 1/k = 9.868 (Set-I) and 
1/k = 8.797 (Set-II) in Table [| correspond approximately to the charm quark. Therefore we 
study the meson correlators with these parameters in the successive calculations. 

Figure ^| shows the hyperfine splitting for charmonium [rrij/^ — m^). For comparison, 
results of other group are shown simultaneously. Here we notice that these results largely 
depend on how to determine the lattice cutoff scale. Results of Fermilab action on the isotropic 
lattice pi]], whose scale is determined from the physical value of the string tension, are roughly 
consistent with our results. Results of heavy relativistic action on anisotropic lattice Jl^] and 
NRQCD[p5ll show the similar tendency. These scales are determined from the Sommer scale 
and the IP-IS splitting respectively. However, all the quenched results are roughly a half 
of experimental results 117(2)MeV[]|]. The differences from the experimental value of the 



hyperfine splitting can be partly explained with the dynamical quark effects [26]. 



4.3 Variational analysis 

We examine the variational analysis for constructing more optimized operators. With this 
analysis we can construct series of operators which have better matching to the states of 
interest. 

Here we explain the principle of variational analysis briefly. Firstly, we suppose that the 
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Figure 4: Hyperfine splitting for charmonium. On this figure results of the other group are 
also shown. These groups|l4], 25, 26 used the Sommer scale, the string tension and IP-IS 



splitting for the determination of the lattice cutoff scale respectively. 



state generated by mesonic operator on the lattice is the linear combination of eigenstates for 
Hamiltonian. Practically we assume that the generated states consist of N linearly indepen- 
dent states. Then we prepare N mesonic operators, which have the same quantum numbers. 
These operators are constructed with different smearing function w(x) in Eq.(^I|). Then we 
get the N x N correlator matrix, C^(t), as the following, 

C%(t) = ^{O$\2,t)oW\0)). (24) 

Because of the symmetric nature of this matrix, we can get the diagonalized correlators which 
are optimized correlators for the state of interest. 

In our study this variational analysis is applied in the minimal space including the ground, 
such as 1st excited and 2nd excited state. We have to prepare the smearing functions so that 
this analysis works well in this condition. Thus we construct the smearing functions 4>i(r) 
using the Schrodinger equation with the potential model, 



1 d 2 1(1 + 1) _ , . 



2mR dr 2 2m#r 2 



yi(r) = Eyi(r), 

yi(r) = r<pi(r), (25) 



where V(r) is the static quark potential measured on our lattice and niR = 1.5/2 GeV. The 
spin interaction is neglected and we calculate only S-state (1 = 0). Figure || shows the lowest 
three solutions (f>o(r) and the measured wave function at T = for Set-II. 

We calculate the diagonalized correlator using the three types of smearing function in 
Fig. [|. The orthogonal matrices are obtained at each t, then we adopt the averaged one as 
the orthogonal matrix which is used for the calculation of the diagonalized correlator. With 
the orthogonal matrix we obtain the diagonalized correlator. Then its effective mass plots 
(Set-II) are shown in Fig. o. We fit the these data at t = 10-36 for Set-I and t = 14-48 for 
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Figure 5: Smearing functions for the variational analysis for Set-II. These are calculated with 
the Schrodinger equation with the potential model. The circle and square symbols are the wave 
function measured on the lattice at T = for the pseudoscalar and the vector respectively. 





Ps(2S) 


V(2S) 


771 Ps(2S) — ^Ps(lS) 


m V(2S) ~ m V(lS) 


Set-I 


0.7688(34) 


0.7794(38) 


0.0868(30) 


0.0848(33) 


(MeV) 


3460(155) 


3507(157) 


391(22) 


382(23) 


Set-II 


0.5641(66) 


0.5696(73) 


0.0913(64) 


0.0869(70) 


(MeV) 


3587(55) 


3622(58) 


581(41) 


553(45) 


Exp. value (MeV) |3] 


3594(5) 


3685.96(9) 


614(5) 


589.07(13) 



Table 6: Spectroscopy for 2S state and 2S-1S splitting of charmonium. For each sets the 
upper values are in a T unit and the lower one are in physical unit. The error in physical unit 
includes the error of a T . 

Set-II, and get the results summarized in Table ||. 

The extracted mass of IS state are consistent with previous results in Table |j| In the 
Table || results are presented in a T unit and physical unit for each sets, where the latter case 
includes the error of a r . 

The results of 2S-1S splitting are consistent with the experimental values within a statis- 
tical error. This is in contrast to the case of the hyperfine splitting in which we obtain the 
half of the experimental value, 2S-1S splitting is in good agreement with experimental value. 

The variational analysis can directly extract one of correlators for the ground and excited 
states. Therefore this analysis is useful to investigate the excited state of hadron. In the next 
section we mainly use this analysis at T > 0. 



5 Results at Finite temperature 
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Figure 6: Effective mass of diagonalized correlator for Set-II at T = 0. The left figure is the 
result of pseudoscalar meson and the right one is that of vector meson. Only ground and 1st 
excited one are shown because of the large noise for the 2nd excited state. 

5.1 Temperature dependence of correlators 

In accordance with our strategy the optimized operators which are constructed at T = are 
applied to the case at T > 0. The optimized operators are defined from the variational analysis 
in the previous section. We construct the correlator matrix at T > with the same basis ( 
smearing function ) as T = ( Fig. ||). Then the orthogonal matrix defined at T = are 
operated to this correlator matrix, and its diagonal correlators are discussed in this subsection. 
Here we call these correlators as Ci(i), C^t) and C^{t) in the order of their effective masses. 
Figure [7| shows the effective masses of C\(t) and C2(t) for Set-II, although the effective mass 
of Cs(t) is too noisy. 

Below T c ( Nt = 26 ) we can find the plateau of effective mass. The ground and the excited 
states seem to be observed below T c . The masses of these states are almost the same as T = 
or slightly larger. Therefore we conclude that these correlators have little thermal effects at 
this temperature, and the spectral structure seems to keep the form at T = 0. However it 
is difficult to identify the plateau precisely and determine the mass quantitatively with the 
present statistics. More detailed analysis with higher statistics may open a stage to discuss 
the potential mass shift of charmonium near to the T c f§. 

Above T c ( Nt = 22 and 16 ) effective masses have no clear plateau in whole t region. 
These behaviors at least signal significant change of correlators when the system crosses T c . 
This behavior appears noticeably in the effective mass of C2(t). In the case of the light quark 
system investigated in the Ref . || , the effective masses increase as T in the pseudoscalar and 
the vector channels. The observed behavior in present work, however, shows qualitatively 
different nature of the correlators. 

As the comparison with above results, Fig. || shows the effective mass plots at T > for 
the correlators with smeared source by the wave function and point sink, described in Sect. 4.1 

These results are consistent with the variational analysis. The remarkable change in the 
vector channel is apparently seen. Then the order of effective masses for the pseudoscalar and 
the vector mesons is reversed above T c , the same as the case of free quark. This phenomenon 
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Figure 7: Effective mass of diagonalized correlator for Set-II at T > 0. Orthogonal matrix 
at T = are used for diagonalization. The left figure is the result of the pseudoscalar meson 
and the right one is that of the vector meson. Because of the large noises only the lowest two 
effective masses are shown. 



is also brought about in the system of light quarks ||. 

From the results of this subsection, we observed that the temporal mesonic correlators 
change drastically when the system goes through the deconfining transition. We can consider 
the several pictures which can explain (consistent with) the change for correlators. One of 
these pictures is that the mesonic bound state disappear in the deconfinement phase. This 
is the most interesting case as stated in our motivations. In order to discuss this picture we 
investigate the correlation between c and c in the next subsection. 



5.2 Wave function 

In this subsection the cc bound state at T > 0, especially in the deconfinement phase, is 
discussed in the light of "wave function" . The definition of the "wave function" in the Coulomb 
gauge is as follows. 

w M (r,t) = Y^{q{x + r,t) lM q{S,t)0^ M {G)) (26) 

X 

Here this definition is the same form as the wave function at T = 0. This wave function shows 
the spatial correlation between q and q, and gives us a hint of the mesonic bound state from 
its t dependence. In the case of free quarks, qq has no bound state, then the wave function 
ought to broaden with t. On the other hand suppose quark and anti-quark form a bound 
state, the wave function holds the stable shape with t. We can discuss the existence of such 
the bound state by observing the i-dependence of the wave function. For this purpose we 
compare the correlation at spatial origin with another spatially separated point at each t. 
Therefore we define the wave function normalized at the spatial origin, (f>M(f,t), as follows, 

<t>M[r,t) = . (27) 

w M {0,t) 

From now on the wave function denotes this normalized definition. 

Since the question is whether wave function has a stable shape or not, it is not necessary 
to use the optimized operator. Therefore the smeared source function with exponential form 
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Figure 8: Effective masses at T > 0. The left figure is that of Set-I and the right one is 
Set-II. In the figures of Set-I, JV t = 20,16 and 12 correspond to T/T c ~ 0.93,1.15 and 1.5 
respectively. In the figures of Set-II, Nt = 26,22 and 16 are T/T c = 0.93, 1.10 and 1.53. The 
symbols with dotted lines show the case of free quarks. 



defined as with Eq . (|23|) , is used for the analysis of i-dependence of the wave function. Fig. |9| 
shows the results at T > with the smeared source function which is slightly wider than the 
observed wave function at T = 0. The wave functions composed of free quark propagators 
are also shown together. 

As is shown in the Fig. ||, the behaviors of the observed wave functions are clearly different 
from that of the free quark case at each temperature and in each mesonic channel. In the 
free quark case the wave functions are broadening as t as expected. On the other hand, the 
observed wave functions are stable with the slightly narrower shape than source function. 
These behaviors are independent of the source function. 

For a visible expression, we define the averaged orbital radius, ro, as 

rl(t) = J23z 2 ^ M (z,t), (28) 

z 

where we suppose a spherical symmetric wave function and the sum is over z axis. These 
ro(t) of Set-I and Set-II are shown in Fig. [u]. Figure 10 shows the t dependence of r$ in the 



physical unit, where the error estimation takes into account the error of a a and a T . These 
results for Set-I and II are roughly consistent with each other. The behavior of observed ro 
are obviously different from free quark case, and shows the stable behavior with respect to t. 
This asymptotic value for the vector channel is larger than that of the pseudoscalar. 

Below T c , ro for each channel are almost same as the case of T = 0. At T ~ 1.5T C the 
observed ro are slightly larger than that of T < T c . However the wave function, even in the 
deconfinement phase, seem to have the stable ro independently of the source function. There- 
fore we conclude the same strong spatial correlation as below T c survives in the deconfinement 
phase at each mesonic channel. 

The results of this subsection seem to suggest the presence of the hadronic state in the 
deconfinement phase, at least up to 1.5T C . This is the opposite situation to the picture 
mentioned in the previous subsection. 
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Figure 9: The results of "wave function" : dotted lines (full symbols) show free quark case, 
and dashed line (no symbol) show the initial shape of "wave function" i.e. shape of source 
function. The different symbols correspond to the different t step. The result of pseudoscalar 
is left one, and that of vector is the other. The top figures are for Set-I and the bottom ones 
are for Set-II. 

6 Conclusion 

In this paper we explored the charmonium correlators in the Euclidean temporal direction 
at T > using the quenched lattice QCD simulation. The high resolution in the temporal 
direction is achieved by employing the anisotropic lattice. We examined the thermal effect on 
the correlators based on the following two quantities: the correlators between the optimized 
operators tuned at T = 0, and the i-dependence of qq wave function. 

In the calculation of the wave function in the Coulomb gauge, we found the strong spatial 
correlation even above T c , up to at least 1.5T C . These results indicate the quark and the 
anti-quark tend to be close each other even in the deconfinement phase. The similar result is 
reported in the light quark mass region in the Ref. ||. In the case of charm quark mass these 
results suggest that the J /if) may not easily be resolved until ~ 1.5T C . On the other hand, 
we also observed significant change of the nature of the correlators between the operators 
optimized at T = 0. This was signaled by the drastic change of the behavior of effective mass. 
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Figure 10: t dependence of r$. Left figure is for Set-I and right one is for Set-II. 

This situation is interesting, and at the same time puzzling. It is possible to consider several 
pictures which are able to explain our results. For example, the mesonic spectral function 
still have some peaks above T c , and its width are broadened with thermal effects. Such a 
situation naturally explains the observed results in this work. The existence of hadronic 
modes just above T c were also suggested by the previous works||. From the analysis of 
temporal meson correlators at T > 0, the QCD vacuum still has non-perturbative nature 
above the phase transition, and far from the perturbative plasma state of quark and gluon. In 
spite of several systematic uncertainties our results in the present work have important and 
interesting implications on the fate of hadronic states above the critical temperature. 

Our approach is directly applicable to the dynamical configuration. If the full QCD 
simulation on the anisotropic lattice is appropriately implemented. Then it is interesting to 
investigate what effect our results receive from the dynamical quarks. To achieve the high 
resolution in the temporal direction, we adopted the anisotropic lattice and employed the O(a) 
improved Wilson quark action on it. These implementation also useful to study the heavy 
particle such as the glueballs |]27[| and scalar mesons. For the correlators of these states we are 
inevitably forced to extract the signals at the short time separation. The large temporal lattice 
cutoff enables us to simulate the heavy particles on lattices of moderate size with keeping the 
finite lattice cutoff effect small. The detailed information in the temporal direction is also 
significantly useful for the direct extraction of the spectral function from the lattice data. This 
approach may give us the further information on the spectral structure of mesons at T > 0. 
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